## PS Replication
## Farris, Key, and Sumner

setwd("C:/Users/Jane Sumner/Dropbox/Research Money Project/Submissions/PS_negotiation/Revision/Materials Resubmitted/Conditional Acceptance/Replication Files")


rep.data <- read.csv("replication-data.csv")
rep.data$highestdegree <- factor(rep.data$highestdegree,
                                levels=c("Bachelor's degree","Master's degree","Doctoral degree"))


## descriptives
table(rep.data$highestdegree)/sum(table(rep.data$highestdegree))
table(rep.data$publicprivate)/sum(table(rep.data$publicprivate))

# percent receiving research funds:
table(rep.data$researchfunds)/sum(table(rep.data$researchfunds))*100

# median research funds:
summary(rep.data$totalmoney2[which(rep.data$totalmoney2>0)])

# median research funds by institution type:
tapply(rep.data$totalmoney2[which(rep.data$totalmoney2>0)],rep.data$inst_type[which(rep.data$totalmoney2>0)],summary)

# type of leave conditional upon receiving one:
table(rep.data$researchleave[which(rep.data$researchleave!="No")])/sum(table(rep.data$researchleave[which(rep.data$researchleave!="No")]))

## SE of proportion
sep <- function(p,N){return(sqrt((p*(1-p))/N))}


### Figure 1
fig1 <- rep.data[,c("moving","housingsearch","summer","tenureclock","spousal","visa","mortgage","defer")]
p <- apply(fig1,2,function(x){
  sum(grepl("Yes",x),na.rm=T)/sum(table(x))
})
N <- apply(fig1,2,
           function(x){
             sum(table(x))
           })
seps <- sep(p,N)
outputs <- data.frame(p,N,seps)
outputs$lb <- outputs$p-1.96*outputs$sep
outputs$ub <- outputs$p+1.96*outputs$sep
outputs$name <- c("Moving",
                  "Housing Search",
                  "Summer 1/9th",
                  "Tenure Clock",
                  "Spousal Hire",
                  "Visa Assistance",
                  "Mortgage Assistance",
                  "Postdoc Deferral")
outputs$name <- gsub(" ","\n",outputs$name)
outputs <- outputs[order(outputs$p),]

par(mfrow=c(3,2))
par(mar=c(3.5,9,2,1))
plot(outputs$p,1:nrow(outputs),pch=16,xlim=c(0,1),xlab="",ylab="",yaxt="n",
     main="Proportion Receiving Compensation and Personal Support Items")
abline(v=seq(0,1,.05),lty=2,col="gray")
points(outputs$p,1:nrow(outputs),pch=16,cex=1.5)
segments(y0=1:nrow(outputs),x0=outputs$lb,x1=outputs$ub,lwd=2)
axis(2,at=1:nrow(outputs),labels=outputs$name,las=2,cex.axis=1)
title(xlab="Proportion of Sample", line=2)

fig1 <- rep.data[,c("computer","software","furniture","space")]
p <- apply(fig1,2,function(x){
  sum(grepl("Yes",x),na.rm=T)/sum(table(x))
})
N <- apply(fig1,2,
           function(x){
             sum(table(x))
           })
seps <- sep(p,N)
outputs <- data.frame(p,N,seps)
outputs$lb <- outputs$p-1.96*outputs$sep
outputs$ub <- outputs$p+1.96*outputs$sep
outputs$name <- c("Computer",
                  "Software",
                  "Furniture",
                  "Lab Space")
outputs$name <- gsub(" ","\n",outputs$name)
outputs <- outputs[order(outputs$p),]

plot(outputs$p,1:nrow(outputs),pch=16,xlim=c(0,1),xlab="",ylab="",yaxt="n",
     main="Proportion Receiving General Support Items ")
abline(v=seq(0,1,.05),lty=2,col="gray")
points(outputs$p,1:nrow(outputs),pch=16,cex=1.5)
segments(y0=1:nrow(outputs),x0=outputs$lb,x1=outputs$ub,lwd=2)
axis(2,at=1:nrow(outputs),labels=outputs$name,las=2,cex.axis=1)
title(xlab="Proportion of Sample", line=2)

fig1 <- rep.data[,c("researchfunds","travel","coursereduction","RAs","researchleave","bookworkshop","speakers","subscriptions")]
p <- apply(fig1,2,function(x){
  sum(grepl("Yes",x),na.rm=T)/sum(table(x))
})
N <- apply(fig1,2,
           function(x){
             sum(table(x))
           })
seps <- sep(p,N)
outputs <- data.frame(p,N,seps)
outputs$lb <- outputs$p-1.96*outputs$sep
outputs$ub <- outputs$p+1.96*outputs$sep

outputs$name <- c("Research Funds",
                  "Travel Funding",
                  "Course Reduction",
                  "RA money",
                  "Research Leave",
                  "Book Workshop",
                  "Conference Host",
                  "Subscription")
outputs$name <- gsub(" ","\n",outputs$name)
outputs <- outputs[order(outputs$p),]

plot(outputs$p,1:nrow(outputs),pch=16,xlim=c(0,1),xlab="",ylab="",yaxt="n",
     main="Proportion Receiving Research Support Items ")
abline(v=seq(0,1,.05),lty=2,col="gray")
points(outputs$p,1:nrow(outputs),pch=16,cex=1.5)
segments(y0=1:nrow(outputs),x0=outputs$lb,x1=outputs$ub,lwd=2)
axis(2,at=1:nrow(outputs),labels=outputs$name,las=2,cex.axis=1)
title(xlab="Proportion of Sample", line=2)

fig1 <- rep.data[,c("books","limitprep","TAs","classpriority")]
p <- apply(fig1,2,function(x){
  sum(grepl("Yes",x),na.rm=T)/sum(table(x))
})
N <- apply(fig1,2,
           function(x){
             sum(table(x))
           })
seps <- sep(p,N)
outputs <- data.frame(p,N,seps)
outputs$lb <- outputs$p-1.96*outputs$sep
outputs$ub <- outputs$p+1.96*outputs$sep
outputs$name <- c("Books",
                  "Limit Prep",
                  "Teaching Assistant",
                  "Class Priority")
outputs$name <- gsub(" ","\n",outputs$name)
outputs <- outputs[order(outputs$p),]

plot(outputs$p,1:nrow(outputs),pch=16,xlim=c(0,1),xlab="",ylab="",yaxt="n",
     main="Proportion Receiving Teaching Support Items ")
abline(v=seq(0,1,.05),lty=2,col="gray")
points(outputs$p,1:nrow(outputs),pch=16,cex=1.5)
segments(y0=1:nrow(outputs),x0=outputs$lb,x1=outputs$ub,lwd=2)
axis(2,at=1:nrow(outputs),labels=outputs$name,las=2,cex.axis=1)
title(xlab="Proportion of Sample", line=2)

fig1 <- rep.data[,c("membership","servicereduc")]
p1 <- sum(grepl("Yes",fig1$membership))/sum(table(fig1$membership))
p2 <- sum(fig1$servicereduc=="Yes, but the exact terms were unspecified.",na.rm=T)/sum(table(fig1$servicereduc))
p3 <- sum(fig1$servicereduc %in% c("Yes, in some other way.","Yes, no service for some or all years beyond the first year.",
                                   "Yes, no service in the first year."),na.rm=T)/sum(table(fig1$servicereduc))
p <- c(p1,p2,p3)
N <- apply(fig1,2,
           function(x){
             sum(table(x))
           })
N <- c(N,N[2])
seps <- sep(p,N)
outputs <- data.frame(p,N,seps)
outputs$lb <- outputs$p-1.96*outputs$sep
outputs$ub <- outputs$p+1.96*outputs$sep
outputs$name <- c("Professional Memberships",
                  "Informal Service Reduction",
                  "Formal Service Reduction")
outputs$name <- gsub(" ","\n",outputs$name)
outputs <- outputs[order(outputs$p),]

plot(outputs$p,1:nrow(outputs),pch=16,xlim=c(0,1),xlab="",ylab="",yaxt="n",
     main="Proportion Receiving Service and Professional Development Items ")
abline(v=seq(0,1,.05),lty=2,col="gray")
points(outputs$p,1:nrow(outputs),pch=16,cex=1.5)
segments(y0=1:nrow(outputs),x0=outputs$lb,x1=outputs$ub,lwd=2)
axis(2,at=1:nrow(outputs),labels=outputs$name,las=2,cex.axis=1)
title(xlab="Proportion of Sample", line=2)


## Figure 2

percentage.tests <- function(variable,group,type,name){
  variable2 <- ifelse(grepl("Yes",variable[which(group==type)]),1,0)
  p <- table(variable2)/sum(table(variable2))
  N <- sum(table(variable2))
  seps <- sep(p,N)
  if(length(p)==2){
    up <- p + 1.96*seps
    low <- p -1.96*seps
  } else{
    p <- 0
    up <- 0
    low <- 0
  }
  return(list(p=data.frame(name=name,p=p[2],lb=low[2],ub=up[2],N=N)))
}

to.plot <- function(type,plotname){
  df <- rbind(percentage.tests(rep.data$anymoney,rep.data$inst_type,type,"Money")$p,
              percentage.tests(rep.data$books,rep.data$inst_type,type,"Books")$p,
              percentage.tests(rep.data$bookworkshop,rep.data$inst_type,type,"Book Workshop")$p,
              percentage.tests(rep.data$classpriority,rep.data$inst_type,type,"Class Priority")$p,
              percentage.tests(rep.data$computer,rep.data$inst_type,type,"Computer")$p,
              percentage.tests(rep.data$conference,rep.data$inst_type,type,"Travel")$p,
              percentage.tests(rep.data$coursereduction,rep.data$inst_type,type,"Course Reduction")$p,
              percentage.tests(rep.data$defer,rep.data$inst_type,type,"Defer for Post-Doc")$p,
              percentage.tests(rep.data$furniture,rep.data$inst_type,type,"Furniture")$p,
              percentage.tests(rep.data$housingsearch,rep.data$inst_type,type,"Housing Search")$p,
              percentage.tests(rep.data$limitprep,rep.data$inst_type,type,"Limit Prep")$p,
              percentage.tests(rep.data$membership,rep.data$inst_type,type,"Memberships")$p,
              percentage.tests(rep.data$mortgage,rep.data$inst_type,type,"Mortgage Assistance")$p,
              percentage.tests(rep.data$moving,rep.data$inst_type,type,"Moving Expenses")$p,
              percentage.tests(rep.data$RAs,rep.data$inst_type,type,"RAs")$p,
              percentage.tests(rep.data$researchleave,rep.data$inst_type,type,"Research Leave")$p,
              percentage.tests(rep.data$servicereduc,rep.data$inst_type,type,"Service Reduction")$p,
              percentage.tests(rep.data$software,rep.data$inst_type,type,"Software")$p,
              percentage.tests(rep.data$space,rep.data$inst_type,type,"Lab Space")$p,
              percentage.tests(rep.data$speakers,rep.data$inst_type,type,"Speakers/Conf. Funds")$p,
              percentage.tests(rep.data$spousal,rep.data$inst_type,type,"Spousal Hire")$p,
              percentage.tests(rep.data$subscriptions,rep.data$inst_type,type,"Subscriptions")$p,
              percentage.tests(rep.data$summer,rep.data$inst_type,type,"Summer 1/9th")$p,
              percentage.tests(rep.data$tenureclock,rep.data$inst_type,type,"Tenure Clock")$p,
              percentage.tests(rep.data$TAs,rep.data$inst_type,type,"TAs")$p,
              percentage.tests(rep.data$visa,rep.data$inst_type,type,"Visa Assistance")$p)
  #omitted <- paste0(df$name[is.na(df$p)],collapse=", ")
  df$p <- ifelse(is.na(df$p),0,df$p)
  
  df$lb <- ifelse(df$lb<0,0,df$lb)
  df$ub <- ifelse(df$ub>1,1,df$ub)
  df <- df[order(df$p),]
  df$order <- order(df$p)
  
  par(mar=c(3.5,8,2,1))
  plot(df$p,df$order,pch=16,xlim=c(0,1),xlab="",ylab="",yaxt="n",main=paste0(plotname," (N = ",unique(df$N),")"))
  abline(v=seq(0,1,.05),lty=2,col="gray")
  points(df$p,df$order,pch=16)
  segments(y0=df$order,x0=df$lb,x1=df$ub)
  axis(2,at=df$order,df$name,las=2,cex.axis=.7)
  title(xlab="Proportion of Sample", line=2)
  
  return(list(df=df))
}

par(mfrow=c(3,2))
to.plot("Bachelor's degree-Public","Bachelor's Degree - Public")
to.plot("Bachelor's degree-Private","Bachelor's Degree - Private")

to.plot("Master's degree-Public","Master's Degree - Public")
to.plot("Master's degree-Private","Master's Degree - Private")
to.plot("Doctoral degree-Public","Doctoral Degree - Public")
to.plot("Doctoral degree-Private","Doctoral Degree - Private")
par(mfrow=c(1,1))

## end of code.
